8 Porous media flow
With the fundamental process model for incompressible flow of a Newtonian fluid, we can now investigate a number of special situations. Our first focus will be on the flow of a Newtonian fluid through a porous background medium, and result in the derivation of Darcy’s law.
8.1 Darcy’s law
We assume a Newtonian fluid flowing through a porous medium. Experiments show the following: The total discharge \(Q\) across an area \(A\) within the porous medium is proportional
- to the size of the cross-sectional area \(A\), and
- to the pressure difference \(p_0 - p_1\),
It is inversely proportional
- to the length of the flow path \(L\), subject to the pressure difference, and
- to the fluid’s viscosity \(\mu\).
The proportionaly constant is denoted by \(\kappa\) and referred to as the intrinsic permeability. The relation hence reads \[ \begin{aligned} Q = \kappa \, A \, (p_0 - p_1) \, \mu^{-1} \,L^{-1} = - \kappa \frac{A (p_1 - p_0)}{\mu L} \end{aligned} \]
Regrouping and considering the limit \(\lim \, L \to 0\) yields \[ \begin{aligned} \underbrace{\frac{Q}{A}}_{=q} = - \frac{\kappa}{\mu} \nabla p, \end{aligned} \]
in which \(q\) is the specific discharge per unit area. How can we understand this?
An ad-hoc argument is to assume the internal friction \(\mu \triangle \mathbf v\) in the small Reynolds Number regime (Stokes flow) to be given by
\[ \begin{aligned} \mu \triangle \mathbf v \approx - \mu \frac{\phi \mathbf v}{\kappa}. \end{aligned} \tag{8.1}\]
The internal stress is hence assumed to be proportional to the porosity \(\phi\), to the velocity \(\mathbf v\) and inversely proportional to the intrinsic permeability \(\kappa\). Substitution yields
\[ \begin{aligned} \underbrace{\phi \mathbf v}_{= \mathbf q} = - \frac{\kappa}{\mu} \left( \nabla p - \rho \mathbf g \right), \label{eq:darcy} \end{aligned} \]
in which we again identify specific discharge \(\mathbf q\).
8.2 Mathematical homogenization
Homogenization is a technique applied to differential equations that describe physical phenomena associated with a domain exhibiting heterogeneities and/or geometric features at two scales or more, typically the macro-scale and the micro-scale. Homogenization allows to consider processes at different scales. The term mathematical theory of homogenization is due to Babuska (1975).
Here, we will use homogenization to derive a continuum mechanical process model for flow through a porous medium referred to as Darcy’s equations.
We assume that the porous medium has some periodic structure on the micro-scale. This means that there exists a micro-scale region \(Y\) of characteristic length \(l\) within the macro- scale domain \(X\) of characteristic length \(L\). The ratio \(\epsilon = l/L\) is hence a small parameter of the system.
The major trick is now to introduce a scaled coordinate
\[y=\frac{x}{\epsilon}\]
that allows us to navigate within the small-scale domain, while having values of reasonable magnitude. Assuming \(L=1\), which is always possible, results in \(\epsilon=l\), and hence \(0 \leq y \leq 1\) in a \(Y\)-cell.
Any function \(f^\epsilon : X \rightarrow \mathbb R^d\) that is affected by processes on the micro-scale can now be interpreted as a two-scale function
\[f : X \times Y \rightarrow \mathbb R^d, (x,y) \rightarrow f(x,y).\]
We furthermore assume that \(f\) can be decomposed into a power series expansion according to
\[ \begin{aligned} f^\epsilon(x) = f(x,y) &= \sum_{i=0}^n \epsilon^i f^{(i)}(x,y) \\ &= f^{(0)}(x,y) + \epsilon^1 f^{(1)}(x,y) + \epsilon^2 f^{(2)}(x,y) + \mathcal O (\epsilon^3) \end{aligned} \]
8.3 An idealized example
Before deriving Darcy’s model, we look at an idealized example that shows the necessary steps for homogenization 1.
Consider a stationary diffusive process, hence a second order differential equation of unknown variable \(u\) and parameter \(a^\epsilon(x)\) that fluctuates on the micro-scale in a periodic manner:
\[ \begin{aligned} \frac{d}{dx} \left( a^\epsilon(x) \frac{d}{dx} u^\epsilon (x) \right) & = 0 \qquad \text{for} \quad 0 \leq x \leq L\\ a^\epsilon(x) &= \frac{1}{1 + 2 \sin^2 (\pi x / \epsilon)} \end{aligned} \tag{8.2}\]
We need two boundary conditions to solve the system:
\[ u^\epsilon (0) = 0 \qquad \text{and} \qquad u^\epsilon (1) = 1 \]
Now, we apply the homogenization trick: \(u^\epsilon (x)\) is interpreted as a two-scale function in \(x\) and \(y=x/\epsilon\) and furthermore written as a power series in \(\epsilon\):
\[ u^\epsilon(x) = u(x,y) = u^{(0)}(x,y) + \epsilon u^{(1)}(x,y) + \epsilon^2 u^{(2)}(x,y) + \mathcal O (\epsilon^3) \]
Differentiation with respect to \(x\) yields
\[ \begin{aligned} \frac{d}{dx} u(x,y) &= \frac{d}{dx} u\left(x,y(x)\right) \\ &= \partial_x u(x,y) + \partial_y u(x,y) \underbrace{\frac{d}{dx}y(x)}_{= \epsilon^{-1}} \\ &= \partial_x u + \epsilon^{-1} \partial_y u \end{aligned} \tag{8.3}\]
This relation holds also for any of the contributing terms \(u^{(i)}(x,y)\) in the power series expansion. By observing that parameter \(a^\epsilon(x)\) depends on \(y\) only
\[ a^\epsilon(x) = \frac{1}{1 + 2 \sin^2 (\pi x / \epsilon)} = \frac{1}{1 + 2 \sin^2 (\pi y)} = a(y), \tag{8.4}\]
we have the necessary building blocks to homogenize the process model. First, we expand according to Equation 8.3. Note, that \(u^{(i)} := u^{(i)}(x,y)\) is introduce to improve readability. This yields \[ \begin{align*} & \quad \frac{d}{dx} \left( a(y) \frac{d}{dx} \left(u^{(0)} + \epsilon u^{(1)} + \epsilon^2 u^{(2)}\right) \right) \\ & = \frac{d}{dx} \left( a(y) \left(\partial_x u^{(0)} + \epsilon^{-1} \partial_y u^{(0)} + \epsilon \partial_x u^{(1)} + \partial_y u^{(1)} + \epsilon^2 \partial_x u^{(2)} + \epsilon \partial_y u^{(2)}\right) \right)\\ % first line & = \underbrace{\partial_x \left( a(y) \, \partial_x u^{(0)} \right)}_{\text{order} \; \epsilon^0} + \underbrace{\epsilon^{-1}\,\partial_y \left( a(y) \, \partial_x u^{(0)} \right)}_{\text{order} \; \epsilon^{-1}} + \underbrace{\epsilon^{-1}\,\partial_x \left( a(y) \, \partial_y u^{(0)} \right)}_{\text{order} \; \epsilon^{-1}} + \underbrace{\epsilon^{-2}\,\partial_y \left( a(y) \, \partial_y u^{(0)} \right)}_{\text{order} \; \epsilon^{-2}} +\\ % second line & + \underbrace{\epsilon \, \partial_x \left( a(y) \, \partial_x u^{(1)} \right)}_{\text{order} \; \epsilon^1} + \underbrace{\partial_y \left( a(y) \, \partial_x u^{(1)} \right)}_{\text{order} \; \epsilon^{0}} + \underbrace{\partial_x \left( a(y) \, \partial_y u^{(1)} \right)}_{\text{order} \; \epsilon^{0}} + \underbrace{\epsilon^{-1}\,\partial_y \left( a(y) \, \partial_y u^{(1)} \right)}_{\text{order} \; \epsilon^{-1}} +\\ % third line & + \underbrace{\epsilon^2 \, \partial_x \left( a(y) \, \partial_x u^{(2)} \right)}_{\text{order} \; \epsilon^2} + \underbrace{\epsilon \,\partial_y \left( a(y) \, \partial_x u^{(2)} \right)}_{\text{order} \; \epsilon^{1}} + \underbrace{\epsilon \,\partial_x \left( a(y) \, \partial_y u^{(2)} \right)}_{\text{order} \; \epsilon^{1}} + \underbrace{\partial_y \left( a(y) \, \partial_y u^{(2)} \right)}_{\text{order} \; \epsilon^{0}} + \end{align*} \]
The model has to be fullfilled on every scale, so it makes sense to arrange terms according to their order in \(\epsilon\):
\[ \begin{align*} \mathcal O \left(\epsilon^{-2}\right): \quad &\partial_y \left( a(y) \, \partial_y u^{(0)} \right) = 0 && \tag*{🟦} \\ \mathcal O \left(\epsilon^{-1}\right): \quad &\partial_x \left( a(y) \, \partial_y u^{(0)} \right) + \partial_y \left( a(y) \, \partial_x u^{(0)} \right) + \partial_y \left( a(y) \, \partial_y u^{(1)} \right) = 0 && \tag*{🟧}\\ \mathcal O \left(\epsilon^{0}\right): \quad &\partial_x \left( a(y) \, \partial_x u^{(0)} \right) + \partial_y \left( a(y) \, \partial_x u^{(1)} \right) + \partial_x \left( a(y) \, \partial_y u^{(1)} \right) + \partial_y \left( a(y) \, \partial_y u^{(2)} \right) = 0 && \tag*{🟨}\\ \mathcal O \left(\epsilon^{1}\right): \quad &\partial_x \left( a(y) \, \partial_x u^{(1)} \right) + \partial_y \left( a(y) \, \partial_x u^{(2)} \right) + \partial_x \left( a(y) \, \partial_y u^{(2)} \right) + \partial_y \left( a(y) \, \partial_y u^{(3)} \right) = 0 && \tag*{🟩}\\ & \vdots \nonumber \end{align*} \tag{8.5}\]
We hence derived a cascade of differential equations Equation 8.5 \(🟦, 🟧, 🟨, 🟩\), which all have to be fulfilled individually.
Again, boundary conditions are needed. The original boundary conditions are assigned to the leading order terms. Higher order terms take zero boundary conditions:
\[ \begin{align*} &u^{(0)} (0,y) = 0 \qquad u^{(0)} (1,y) = 1 \\ &u^{(i)} (0,y) = 0 \qquad u^{(i)} (1,y) = 0 \qquad &\text{for all} \; i \geq 1 \end{align*} \]
Periodicity on the micro-scale requires that
\[ \begin{align*} &u^{(i)} (x,0) = u^{(i)} (x,1) \qquad &\text{for all} \; i \geq 0 \tag*{🟦}\\ &\partial_y u^{(i)} (0,y) = \partial_y u^{(i)} (1,y) = 0 \qquad &\text{for all} \; i \geq 0 \tag*{🟧} \end{align*} \tag{8.6}\]
Now step by step, the cascade of equations Equation 8.5 \(🟦, 🟧, 🟨, 🟩\) can be solved:
Step 1: We observe that Equation 8.5 \(🟦\) holds, if \(u^{(0)}\) is independent of \(y\).
\[ u^{(0)}(x,y) = u^{(0)}(x) \tag{8.7}\]
Step 2: Substituting the latter into Equation 8.5 \(🟧\) yields
\[ \begin{align*} \underbrace{\partial_y \left( a(y) \, \partial_x u^{(0)} \right)}_{\left( \partial_y a(y) \right) \partial_x u^{(0)} + \underbrace{a(y) \partial_{yx} u^{(0)}}_{=0}} + \partial_y \left( a(y) \, \partial_y u^{(1)} \right) = 0. \end{align*} \tag{8.8}\]
With \(u^{(1)}\) defined according to \[ \begin{align*} u^{(1)}(x,y) = \left( \frac{\int_0^y \frac{1}{a(\tilde y)} d\tilde y}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} - y + c\right) \partial_x u^{(0)} + f(x) \end{align*} \tag{8.9}\]
equation Equation 8.8 is fulfilled.
Step 3: Integrate Equation 8.5 \(🟨\) over the periodic cell \(Y\):
\(\int_0^1\) Equation 8.5 \(🟨\) \(d y\) = \[ \begin{align*} & \int_0^1 \partial_x \left( a(y) \partial_x u^{(0)} \right) d y + \underbrace{\int_0^1 \partial_y \left( a(y) \, \partial_x u^{(1)} \right) d y}_{[a(y) \partial_x u^{(1)}]_0^1 = 0} \\ & + \int_0^1 \partial_x \left( a(y) \partial_y u^{(1)} \right) d y + \underbrace{\int_0^1 \partial_y \left( a(y) \partial_y u^{(2)} \right) d y}_{[a(y) \partial_y u^{(2)}]_0^1 = 0} \end{align*} \]
where we observe several terms to cancel out due to the micro-scale periodicity Equation 8.6 \(🟦\) and \(🟧\). The following two terms remain:
\(\int_0^1\) Equation 8.5 \(🟨\) \(d y\) = \[ \begin{align*} \int_0^1 \partial_x \left( a(y) \partial_x u^{(0)} \right) + \partial_x \left( a(y) \partial_y u^{(1)} \right) d y \end{align*} \]
Substituting in \(u^{(1)}\) given by Equation 8.9 finally yields
\(\int_0^1\) Equation 8.5 \(🟨\) \(d y\) =
\[ \begin{align*} & \int_0^1 \partial_x \left( a(y) \partial_x u^{(0)} \right) + \partial_x \left( a(y) \left( \left( \frac{\frac{1}{a( y)}}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} - 1 \right) \partial_x u^{(0)} \right) \right) d y \\ = & \int_0^1 \partial_x \left( \frac{1}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} \partial_x u^{(0)} \right) d y = \partial_x \left(\underbrace{\int_0^1 \frac{1}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} d y }_{=:a^{(0)}}\right) \partial_x u^{(0)} \\ = & \frac{d}{dx} \left( a^{(0)} \frac{d}{dx} u^{(0)} \right). \end{align*} \]
This is a rather remarkable results, as it allows to determine the macro-scale variable \(u^{(0)}\) based on a less complex (homogenized) model with a lumped or effective coefficient \(a^{(0)}\) instead of the hard to handle fluctuating \(a^\epsilon (x)\). For our particular example, we get
\[ \begin{align*} \int_0^1 \frac{1}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} d y = \frac{1}{2} \end{align*} \]
and hence
\[ \begin{align*} u^{(0)} = x \qquad \text{and} \qquad u^{(1)} = - \frac{1}{ 4 \pi} \sin \left( 2 \pi y \right). \end{align*} \]
The complete solution reads \[ \begin{align*} u(x,y) = u^{(0)} + \epsilon u^{(1)} = x - \frac{\epsilon}{ 4 \pi} \sin \left( 2 \pi y \right) \end{align*} \]
and we indeed see that \[ \begin{align*} \underset{\epsilon \to 0}{\lim} \, u(x,y) = x, \end{align*} \]
as we initially postulated!
8.4 Derivation of Darcy’s model
Now, we do the same to upscale incompressible Newtonian flow in a porous medium to yield Darcy’s model. It is important, though, to specify the assumptions that this will be based upon:
- The porous medium is non-deformable and stationary.
- The fluid is an incompressible Newtonian fluid of constant density \(\rho\).
- The body force is due to gravitational acceleration only.
- The flow field is at steady-state.
A schematic of the situation has been discussed during class and introduced the following variables:
- \(\Omega\): macro-scale domain, which consists of \(\Omega_l\) (fluid domain) and \(\Omega_s\) (porous medium)
- \(Y\): micro-scale domain, which consists of \(Y_l\) (fluid domain) and \(Y_s\) (porous medium)
- \(\Gamma_e\): fluid part of the \(Y\) boundary
- \(\Gamma_i\): fluid-solid boundary within \(Y\)
- \(\Gamma\): union of all fluid-solid boundaries in \(\Omega\)
- \(\mathbf x = (x_1,x_2,x_3)\): the macro-scale coordinate
- \(\mathbf y = (y_1,y_2,y_3)\): the micro-scale coordinate
- \(L\): characteristic length scale of \(\Omega\)
- \(l\): characteristic length scale of \(Y\)
- \(\epsilon = l/L\)
The void space is fully occupied by the fluid, for which we assume the steady-state incompressible Navier-Stokes equations. The superscript \(\epsilon\) denotes that the respective variable is affected by the micro-scale. For all \(x \in \Omega_l\), we get
\[ \begin{align*} \nabla \cdot \rho \mathbf v^\epsilon (\mathbf x) & = 0 \tag*{🟦}\\ \rho \left( \mathbf v^\epsilon(\mathbf x) \cdot \nabla \right) \mathbf v^\epsilon(\mathbf x) &= - \frac{1}{\rho} \nabla p^\epsilon(\mathbf x) + \mu \triangle \mathbf v^\epsilon(\mathbf x) + \rho \mathbf g. \tag*{🟧} \end{align*} \tag{8.10}\]
No-slip boundary’s are assumed, hence \(\mathbf v^\epsilon = 0\) for \(x \in \Gamma_i\). We introduce
\[\mathbf y := \mathbf x/\epsilon\]
and assume that the dependent variables can be expanded into an asymptotic expansion in \(\epsilon\):
\[ \begin{align*} p^\epsilon(\mathbf x) & = p(\mathbf x,\mathbf y) = p^{(0)}(\mathbf x,\mathbf y) + \epsilon p^{(1)}(\mathbf x,\mathbf y) + \epsilon^2 p^{(2)}(\mathbf x,\mathbf y) + ...\\ \mathbf v^\epsilon(\mathbf x) & = \mathbf v(\mathbf x,\mathbf y) = \epsilon^2 \mathbf v^{(0)}(\mathbf x,\mathbf y) + \epsilon^3 \mathbf v^{(1)}(\mathbf x,\mathbf y) + \epsilon^4 \mathbf v^{(2)}(\mathbf x,\mathbf y) + ... \end{align*} \]
Expansion of Equation 8.10 \(🟦\) and \(🟧\) yields:
Equation 8.10 \(🟦\): \[ \begin{align*} 0 &= \nabla \cdot \left( \rho \left( \epsilon^2 \mathbf v^{(0)} + \epsilon^3 \mathbf v^{(1)} \right) \right) \\ &=\epsilon^2 \nabla_x \cdot \left( \rho \mathbf v^{(0)} \right) + \epsilon^3 \nabla_x \cdot \left( \rho \mathbf v^{(1)} \right) + \epsilon \nabla_y \cdot \left( \rho \mathbf v^{(0)} \right) + \epsilon^2 \nabla_y \cdot \left( \rho \mathbf v^{(1)} \right)\\ \end{align*} \]
Equation 8.10 \(🟧\)-LHS: \[ \begin{align*} & \rho \left( \epsilon^2 \mathbf v^{(0)} + \epsilon^3 \mathbf v^{(1)} \right) \cdot \nabla \left( \epsilon^2 \mathbf v^{(0)} + \epsilon^3 \mathbf v^{(1)}\right) = \nonumber \\ &\epsilon^4 \rho \mathbf v^{(0)} \cdot \nabla_x \mathbf v^{(0)} + \epsilon^5 \rho \mathbf v^{(1)} \cdot \nabla_x \mathbf v^{(0)} + \epsilon^5 \rho \mathbf v^{(0)} \cdot \nabla_x \mathbf v^{(1)}+ \epsilon^6 \rho \mathbf v^{(1)} \cdot \nabla_x \mathbf v^{(1)} \nonumber\\ &+\epsilon^3 \rho \mathbf v^{(0)} \cdot \nabla_y \mathbf v^{(0)} + \epsilon^4 \rho \mathbf v^{(1)} \cdot \nabla_y \mathbf v^{(0)} + \epsilon^4 \rho \mathbf v^{(0)} \cdot \nabla_y \mathbf v^{(1)}+ \epsilon^5 \rho \mathbf v^{(1)} \cdot \nabla_y \mathbf v^{(1)}\\ \end{align*} \]
Equation 8.10 \(🟧\)-RHS: \[ \begin{align*} & - \nabla \left( p^{(0)} + \epsilon p^{(1)} + \epsilon^2 p^{(2)}\right) + \mu \nabla \cdot \nabla \left( \epsilon^2 \mathbf v^{(0)} + \epsilon^3 \mathbf v^{(1)} \right) + \rho \mathbf g = \nonumber \\ & - \nabla_x p^{(0)} - \epsilon \nabla_x p^{(1)} - \epsilon^2 \nabla_x p^{(2)} - \epsilon^{-1}\nabla_y p^{(0)} - \nabla_y p^{(1)} - \epsilon \nabla_y p^{(2)} \nonumber \\ & + \epsilon^2 \mu \nabla_x \cdot \nabla_x \mathbf v^{(0)} + \epsilon^3 \mu \nabla_x \cdot \nabla_x \mathbf v^{(1)} + \epsilon \mu \nabla_x \cdot \nabla_y \mathbf v^{(0)} + \epsilon^2 \mu \nabla_x \cdot \nabla_y \mathbf v^{(1)} \nonumber \\ & + \epsilon \mu \nabla_y \cdot \nabla_x \mathbf v^{(0)} + \epsilon^2 \mu \nabla_y \cdot \nabla_x \mathbf v^{(1)} + \mu \nabla_y \cdot \nabla_y \mathbf v^{(0)} + \epsilon \mu \nabla_y \cdot \nabla_y \mathbf v^{(1)} \end{align*} \]
The equations have to hold on every order, hence every power of \(\epsilon\):
\[ \begin{align*} \mathcal O \left(\epsilon^{-1}\right): \quad & 0 = 0 \label{dem1a} \\ & \nabla_y p^{(0)} = 0 \label{dem1b} \tag*{🟦}\\ \mathcal O \left(\epsilon^{0}\right): \quad & 0 = 0 \label{de0a} \tag*{🟧}\\ & \mu \nabla_y \cdot \nabla_y \mathbf v^{(0)} = \nabla_x p^{(0)} + \nabla_y p^{(1)} - \rho \mathbf g \label{de0b}\tag*{🟨}\\ \mathcal O \left(\epsilon^{1}\right): \quad & \nabla_y \cdot \left( \rho \mathbf v^{(0)} \right) = 0 \label{de1a}\tag*{🟩}\\ & \mu \left( \nabla_y \cdot \nabla_y \mathbf v^{(1)} + \nabla_x \cdot \nabla_y \mathbf v^{(0)} + \nabla_y \cdot \nabla_x \mathbf v^{(0)}\right) = \nabla_x p^{(1)} + \nabla_y p^{(2)} \label{de1b}\tag*{🟪}\\ \mathcal O \left(\epsilon^{2}\right): \quad & \nabla_x \cdot \left( \rho \mathbf v^{(0)} \right) + \nabla_y \cdot \left( \rho \mathbf v^{(1)} \right) = 0 \label{de2a}\tag*{🟫}\\ & ... \end{align*} \tag{8.11}\]
with boundary conditions \(\mathbf v^{(0)} = \mathbf v^{(i)} = 0\) for all \(i\) and \(\mathbf x \in \Gamma, \mathbf y \in \Gamma_i\).
Step by step, we solve the system:
Step 1: \(p^{(0)}\) independent of \(y\) is an admissible solution to Equation 8.11 \(🟦\).
Step 2: Since \(\rho\) is constant, we deduce from Equation 8.11 \(🟩\) that \[ \begin{align*} \nabla_y \cdot \mathbf v^{(0)} &= 0 \end{align*} \]
and from Equation 8.11 \(🟫\) that
\[ \begin{align*} \nabla_x \cdot \mathbf v^{(0)} + \nabla_y \cdot \mathbf v^{(1)} &= 0 \label{stepbsimple} \end{align*} \tag{8.12}\]
Step 3: We define so-called homogenized velocity \(q^{(0)}\):
\[ \begin{align*} \mathbf q^{(0)} = \frac{1}{|Y|} \int_{Y_l} \mathbf v^{(0)}(\mathbf x,\mathbf y) dy \end{align*} \]
with \(|Y|\) denoting the volume of \(Y\).
Integration of Equation 8.12 over the void domain of the periodic \(Y\) cell yields
\(\frac{1}{|Y|} \int_{Y_l}\) Equation 8.12 \(dy\) = \[ \begin{align*} & \underbrace{\frac{1}{|Y|} \int_{Y_l} \nabla_x \cdot \mathbf v^{(0)} dy}_{\frac{1}{|Y|} \nabla_x \cdot \int_{Y_l} \mathbf v^{(0)} dy} + \underbrace{\frac{1}{|Y|} \int_{Y_l} \nabla_y \cdot \mathbf v^{(1)} dy}_{\frac{1}{|Y|} \int_{\partial Y_l} \mathbf v^{(1)} \cdot \mathbf n d\sigma} =\\ & \frac{1}{|Y|} \nabla_x \cdot \mathbf q^{(0)} + \frac{1}{|Y|} \left( \underbrace{\int_{\Gamma_i } \mathbf v^{(1)} \cdot \mathbf n d\sigma}_{=0 \; \text{no-slip}} + \underbrace{\int_{\Gamma_l } \mathbf v^{(1)} \cdot \mathbf n d\sigma}_{=0 \; \text{periodicity}} \right) = 0 \label{homvelzero} \end{align*} \tag{8.13}\]
Step 4: Now, we investigate Equation 8.11 \(🟨\) in greater detail. Let’s make the following Ansatz for \(\mathbf v^{(0)}\)
\[ \begin{align*} \label{decompositionofv} \mathbf v^{(0)} = -\frac{1}{\mu} \mathbf w (\mathbf y) \left( \nabla_x p^{(0)} - \rho \mathbf g \right), \end{align*} \tag{8.14}\]
in which \(\mathbf v^{(0)}, \nabla_x p^{(0)}\) and \(\mathbf g\) are vectors and \(\mathbf w\) is a second order tensor, which has to fulfill \[ \begin{align*} \label{bufferw} \nabla_y \cdot \nabla_y \mathbf w &= \nabla_y \mathbf s - \mathbf I & \mathbf y \in Y_l\\ \nabla_y \cdot \mathbf w &= 0 & \mathbf y \in Y_l\\ \mathbf w &= 0 & \mathbf y \in \Gamma_i\\ \end{align*} \tag{8.15}\]
Step 5: Finally, we integrate Equation 8.14 over \(Y_l\): \[ \begin{align*} \mathbf q^{(0)} &= \frac{1}{|Y|} \int_{Y_l} \left( -\frac{1}{\mu} \mathbf w (\mathbf y) \left( \nabla_x p^{(0)} - \rho \mathbf g \right) \right) dy \\ & = -\frac{1}{\mu} \underbrace{\left( \frac{1}{|Y|} \int_{Y_l} \mathbf w (\mathbf y) dy \right)}_{=: \mathbf \kappa} \left( \nabla_x p^{(0)} - \rho \mathbf g \right) \\ & = - \frac{\mathbf \kappa}{\mu} \left( \nabla_x p^{(0)} - \rho \mathbf g \right) \end{align*} \]
Often, Darcy’s model it is also written in terms of piezometric head
\[h(x):=\frac{p}{\rho g} + z\]
and the hydraulic conductivity
\[\mathbf K := \frac{\rho g}{\mu} \mathbf \kappa.\]
It then has the form
\[ \begin{aligned} \mathbf q & = - \mathbf K \, \nabla h \end{aligned} \tag{8.16}\]
Note, that we used \(\mathbf q\) to abbreviate \(\mathbf q^{(0)}\). Note also, that \(\mathbf K\) is no longer a function of the geometry alone.
Combination with the continuity quation yields
\[ \begin{align*} \label{darcyslawinh} \nabla \cdot \left(\mathbf K \, \nabla h \right) = 0 \end{align*} \tag{8.17}\]
8.5 Effective hydraulic conductivity
We previsously saw that the homogenized continuity equation is given by
\[ \begin{align*} & \nabla \cdot \mathbf q (\mathbf x) = 0 \qquad \text{with} \qquad \mathbf{q}(\mathbf x) := \frac{1}{|Y|} \int_{Y_l} \mathbf v^{(0)}(\mathbf x,\mathbf y), \end{align*} \]
in which \(\mathbf x\) denotes the macro-scale coordinate, \(\mathbf y\), the micro-scale coordinate, and \(\mathbf v^{(0)}\) the leading order term in the velocity expansion.
The macro-scale dynamics is accordingly governed by the following system
\[ \begin{align*} \nabla \cdot \mathbf q (\mathbf x) &= 0 \\ \nabla \cdot \left( \mathbf K (\mathbf x ) \cdot \nabla h (\mathbf x) \right) &=0, \end{align*} \tag{8.18}\]
with \(h(\mathbf x) = h_B ( \mathbf x )\) for \(\mathbf x \in \partial \Omega\).
Let’s assume, we are at field scale, and the porous medium is given as a heterogeneous mosaic of smaller porous media patches, e.g. clay lenses in a rock, or a layered acquifer.
Idea: We want to use use the homogenization technique in order to study the effective flow dynamics at field scale.
Assume that \(\mathbf K\) is goverened by fluctuations on the macro-scale. To make this fact explicit, we write it as \(\mathbf K^\epsilon ( \mathbf x )\) and repeat the homogenization procedure as previously introduced (definition of \(\mathbf y = \mathbf x/\epsilon\) and two-scale interpretaion \(\mathbf K^\epsilon ( \mathbf x ) = \mathbf K ( \mathbf x, \mathbf y )\)). note, however, that now \(\mathbf y\) refers to the macro-scale, while \(\mathbf x\) refers to the mega-scale. Again, \(\mathbf x\) denotes the larger spatial scale of the system.
Now, we write \(h ( \mathbf x )\) as an asymptotic expansion in \(\epsilon\):
\[ \begin{align*} h^\epsilon(\mathbf x) & = h(\mathbf x,\mathbf y) = h^{(0)}(\mathbf x,\mathbf y) + \epsilon h^{(1)}(\mathbf x,\mathbf y) + \epsilon^2 h^{(2)}(\mathbf x,\mathbf y) + ... \end{align*} \]
Substitution into Equation 8.18 yields
\[ \begin{align*} &\nabla \cdot \left( \mathbf K^\epsilon (\mathbf x ) \nabla h (\mathbf x) \right) = \nabla \cdot \left( \mathbf K (\mathbf x ) \nabla \left( h^{(0)} + \epsilon h^{(1)} + \epsilon^2 h^{(2)} \right) \right) =\\ &\underbrace{\nabla_x \cdot \left( \mathbf K \nabla_x h^{(0)} \right)}_{\epsilon^0} + \underbrace{\epsilon \nabla_x \cdot \left( \mathbf K \nabla_x h^{(1)} \right)}_{\epsilon^1} + \underbrace{\epsilon^2 \nabla_x \cdot \left( \mathbf K \nabla_x h^{(2)} \right)}_{\epsilon^1} + \\ &\underbrace{\epsilon^{-1}\nabla_x \cdot \left( \mathbf K \nabla_y h^{(0)} \right)}_{\epsilon^{-1}} + \underbrace{\nabla_x \cdot \left( \mathbf K \nabla_y h^{(1)} \right)}_{\epsilon^0} + \underbrace{\epsilon \nabla_x \cdot \left( \mathbf K \nabla_y h^{(2)} \right)}_{\epsilon^1} + \\ &\underbrace{\epsilon^{-1}\nabla_y \cdot \left( \mathbf K \nabla_x h^{(0)} \right)}_{\epsilon^{-1}} + \underbrace{\nabla_y \cdot \left( \mathbf K \nabla_x h^{(1)} \right)}_{\epsilon^0} + \underbrace{\epsilon \nabla_y \cdot \left( \mathbf K \nabla_x h^{(2)} \right)}_{\epsilon^1} + \\ &\underbrace{\epsilon^{-2}\nabla_y \cdot \left( \mathbf K \nabla_y h^{(0)} \right)}_{\epsilon^{-2}} + \underbrace{\epsilon^{-1}\nabla_y \cdot \left( \mathbf K \nabla_y h^{(1)} \right)}_{\epsilon^{-1}} + \underbrace{\nabla_y \cdot \left( \mathbf K \nabla_y h^{(2)} \right)}_{\epsilon^0} = 0, \end{align*} \]
which can be organized in powers of \(\epsilon\):
\[ \begin{align*} \mathcal O \left(\epsilon^{-2}\right): \quad &\nabla_y \cdot \left( \mathbf K \nabla_y h^{(0)} \right) = 0 \label{effhm2} \tag*{🟦}\\ \mathcal O \left(\epsilon^{-1}\right): \quad &\nabla_x \cdot \left( \mathbf K \nabla_y h^{(0)} \right) + \nabla_y \cdot \left( \mathbf K \nabla_x h^{(0)} \right) + \nabla_y \cdot \left( \mathbf K \nabla_y h^{(1)} \right) = 0 \label{effhm1} \tag*{🟧}\\ \mathcal O \left(\epsilon^{0}\right): \quad &\nabla_x \cdot \left( \mathbf K \nabla_x h^{(0)} \right) + \nabla_x \cdot \left( \mathbf K \nabla_y h^{(1)} \right) + \nabla_y \cdot \left( \mathbf K \nabla_x h^{(1)} \right) + \nabla_y \cdot \left( \mathbf K \nabla_y h^{(2)} \right)= 0 \label{effhm0} \tag*{🟨} \end{align*} \tag{8.19}\]
We solve the system according to the following steps:
Step 1: We observe that \[ \begin{align*} \label{hindy} h^{(0)}(\mathbf x, \mathbf y)= h^{(0)}(\mathbf x) \end{align*} \tag{8.20}\] is a solution to Equation 8.19 \(🟦\).
Step 2: Substitution into Equation 8.19 \(🟧\) yields \[ \begin{align*} \label{effhm1subs} \nabla_y \cdot \left( \mathbf K \nabla_y h^{(1)} \right) = -\nabla_y \cdot \left( \mathbf K \nabla_x h^{(0)} \right) \end{align*} \tag{8.21}\]
Step 3: We make an Ansatz for \(h^{(1)}\), which reads \[ \begin{align*} \label{decomph1} h^{(1)} (\mathbf x, \mathbf y ) = \mathbf r (\mathbf x, \mathbf y ) \cdot \nabla_x h^{(0)} (\mathbf x) + f(\mathbf x), \end{align*} \tag{8.22}\] with a vector \(\mathbf r (\mathbf x, \mathbf y )\) that is periodic in \(Y\), and a scalar function \(f(\mathbf x)\). This indeed simplifies the situation, because rather than solving the system Equation 8.21, it suffices to solve \[ \begin{align*} \label{effhm1subssimpl} \nabla_y \cdot \left( \mathbf K \nabla_y \mathbf r \right) = -\nabla_y \mathbf K \end{align*} \tag{8.23}\]
Step 4: Now, we integrate Equation 8.19 \(🟨\) over \(Y\): \[ \begin{align*} \frac{1}{Y} \int_Y \nabla_x \cdot \left( \mathbf K \nabla_x h^{(0)} \right) \, dy \, + \frac{1}{Y} \int_Y \nabla_x \cdot \left( \mathbf K \nabla_y h^{(1)} \right) \, dy \, + \\ \frac{1}{Y} \int_Y \nabla_y \cdot \left( \mathbf K \nabla_x h^{(1)} \right) \, dy \, + \frac{1}{Y} \int_Y \nabla_y \cdot \left( \mathbf K \nabla_y h^{(2)} \right) \, = \, 0. \end{align*} \] By applying divergence theorem and substituting in Equation 8.22, we get \[ \begin{align*} \frac{1}{Y} \int_Y \nabla_x \cdot \left( \mathbf K \nabla_x h^{(0)} \right) \, dy + \frac{1}{Y} \int_Y \nabla_x \cdot \left( \mathbf K \nabla_y (\mathbf r \cdot \nabla_x h^{(0)}+f) \right) \, dy \, +\\ \underbrace{\int_{\partial Y} \mathbf n_y \cdot \left( \mathbf K \nabla_x h^{(1)} + \mathbf K \nabla_y h^{(2)} \right) \, d\sigma}_{=\, 0 \; \text{due to periodicity condition}} \, = \, 0 \end{align*} \] Now, we pull the macro-scale gradient out of the \(\mathbf y\)-integral, which results in \[ \begin{align*} \label{homeffh} 0 &=\nabla_x \cdot \left( \underbrace{\left( \frac{1}{Y} \int_Y \mathbf K \, dy + \frac{1}{Y} \int_Y \mathbf K \nabla_y \mathbf r \right)}_{=:\mathbf K^{\text{eff}} (\mathbf x)} \nabla_x h^{(0)} \right) \\ & = \nabla_x \cdot \left( \mathbf K^{\text{eff}} \nabla_x h^{(0)} \right) \end{align*} \tag{8.24}\] The latter equations Equation 8.24 is the homogenized version of Darcy’s law at field scale. \(\mathbf K^{\text{eff}}\) is referred to as the effective hydraulic conductivity.
8.6 A layered aquifer
We examine the results of the latter section more closely for a vertically layered aquifer. Layers consists of a homogeneous isotropic material. We assume that \(\mathbf K\) has no large scale trend, hence \(\mathbf K (\mathbf x, \mathbf y) = \mathbf K (\mathbf y) = K(y_i)\). From the previous section we saw that it is essential to determine \(\mathbf r\), which was part of the Ansatz we made for the Darcy velocity.
We write Equation 8.23 componentwise:
\[ \begin{align*} & \partial_{y_1} \left( K(y_1) \, \partial_{y_1} r_1 \right) + \partial_{y_2} \left( K(y_1) \, \partial_{y_2} r_1 \right) + \partial_{y_3} \left( K(y_1) \, \partial_{y_3} r_1 \right) = - \partial_{y_1} K(y_1) \label{caseeq1} \tag*{🟦} \\ & \partial_{y_1} \left( K(y_1) \, \partial_{y_1} r_2 \right) + \partial_{y_2} \left( K(y_1) \, \partial_{y_2} r_2 \right) + \partial_{y_3} \left( K(y_1) \, \partial_{y_3} r_2 \right) = - \underbrace{\partial_{y_2} K(y_1)}_{=0} \label{caseeq2} \tag*{🟧}\\ & \partial_{y_1} \left( K(y_1) \, \partial_{y_1} r_3 \right) + \partial_{y_2} \left( K(y_1) \, \partial_{y_2} r_3 \right) + \partial_{y_3} \left( K(y_1) \, \partial_{y_3} r_3 \right) = - \underbrace{\partial_{y_3} K(y_1)}_{=0}\label{caseeq3} \tag*{🟨} \end{align*} \tag{8.25}\]
The only admissible periodic solutions \(r_i(0) = r_l (0)\) are characterized by \(r_2 = r_3=0\) and \(r_1\) only dependent on \(y_1\).
This reduces Equation 8.25 \(🟦\), \(🟧\), \(🟨\) to
\[ \begin{align*} \partial_{y_1} \left( K(y_1) \, \partial_{y_1} r_1(y_1) \right) &= - \partial_{y_1} K(y_1) \\ \end{align*} \]
from which we deduce the solution
\[ \begin{align*} r_1 (y_1) &= l \frac{\int_0^{y_1} \frac{1}{K(y_1)} d y_1}{\int_0^{l} \frac{1}{K(y_1)} d y_1} - y_1. \end{align*} \tag{8.26}\]
Substitution into Equation 8.24 results in
\[ \begin{align*} K_{11}^{\text{eff}} &= \frac{1}{l} \int_0^l K(y_1) \, d y_1 + \frac{1}{l} \int_0^l K(y_1) \partial_{y_1} r_1 (y_1) \, d y_1 =\frac{l}{\int_0^{l} \frac{1}{K(y_1)} d y_1} \end{align*} \]
and
\[ \begin{align*} K_{22}^{\text{eff}} = K_{33}^{\text{eff}} &= \frac{1}{l} \int_0^l K(y_1) \, d y_1 + \underbrace{\frac{1}{l} \int_0^l K(y_1) \partial_{y_{2/3}} r_{2/3} (y_1) \, d y_1}_{=0}. \end{align*} \]
These results can be combined into a complete effective hydraulic conductivity tensor \[ \begin{align*} \mathbf K^{\text{eff}} = \left(\begin{matrix} \frac{l}{\int_0^{l} \frac{1}{K(y_1)} d y_1} & 0 & 0 \\ 0 & \frac{1}{l} \int_0^{l} K(y_1) d y_1 & 0 \\ 0 & 0 & \frac{1}{l} \int_0^{l} K(y_1) d y_1 \end{matrix} \right) \end{align*} \]
The example is inspired by an introductory chapter in the book of Jacob Bear ↩︎